data = readtable('../../data/Netherlands/Holland_PF_v2.xlsx');
date = data.Var1;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T = length(date(startdate:enddate));

option.detrend = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

time = data.Var1(startdate:enddate);
pi = data.pi(startdate:enddate);
rpcgdpgr = data.x(startdate:enddate);
tau = data.tau(startdate:enddate);
g = data.g(startdate:enddate);
ynom10 = data.y10(startdate:enddate);
surplus = tau - g;

debt_book = data.debt_gdp(startdate:enddate); % book value of debt/gdp ratio
debt = data.debt_gdp_mkt(startdate:enddate);

%%convenience yield from uk
cy = 0.015;
ynom10 = log(1 + ynom10 + cy);

tau = tau + debt * cy;

deltalogg = [log(data.g(startdate)) - log(data.g(startdate - 1)); log(g(2:end)) - log(g(1:end - 1))];
deltalogtau = [log(data.tau(startdate)) - log(data.tau(startdate - 1)); log(tau(2:end)) - log(tau(1:end - 1))];

dy = data.dy(startdate:enddate);
pdm = log(1 ./ dy); % log price/dividend ratio
divgrm = data.divgrm(startdate:enddate); % log dividend growth

x0 = mean(rpcgdpgr);
pi0 = mean(pi);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
y0nom_10 = mean(ynom10);
surplus0 = mean(surplus);

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if option.detrend == 1
    tts(1) = tau(1);
    gts(1) = g(1);

    for t = 2:T
        gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
        tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define the ordering in the VAR
inflpos = 1;
ypos = 2;
gdppos = 3;
divgrpos = 4;
coint_div = 5;
pdpos = 6;
dtpos = 7;
coint_tax = 8;
dgpos = 9;
coint_spending = 10;

%% raw data
clear Xraw;
Xraw(:, inflpos) = pi;
Xraw(:, gdppos) = rpcgdpgr;
Xraw(:, ypos) = ynom10;
Xraw(:, dtpos) = deltalogtau;
Xraw(:, coint_tax) = log(tau);
Xraw(:, dgpos) = deltalogg;
Xraw(:, coint_spending) = log(g);
Xraw(:, pdpos) = pdm;
Xraw(:, divgrpos) = divgrm;
Xraw(:, coint_div) = cumsum(divgrm);

% VAR elements
infl = pi - pi0;
x = rpcgdpgr - x0;
yield10 = ynom10 - y0nom_10;
dg = deltalogg - g0;
dt = deltalogtau - tau0;
pd = pdm - A0m;
% divgr = divgrm - mu_m;

X2(:, inflpos) = infl;
X2(:, ypos) = yield10;
X2(:, gdppos) = x;
X2(:, dtpos) = dt;
X2(:, coint_tax) = log(tau) - mean(log(tau));
X2(:, dgpos) = dg;
X2(:, coint_spending) = log(g) - mean(log(g));
X2(:, divgrpos) = divgrm - mean(divgrm);
X2(:, coint_div) = cumsum(divgrm) - mean(cumsum(divgrm));
X2(:, pdpos) = pd;

% trend adjustment
if option.detrend == 1
    X2(:, coint_tax) = log(tts) - mean(log(tts));
    X2(:, coint_spending) = log(gts) - mean(log(gts));
end

z_var = X2(2:end, :);
z_varlag = X2(1:end - 1, :);
Y = z_var;
X = z_varlag;

N = cols(X2);
I = eye(N);
I_pi = I(:, inflpos);
I_gdp = I(:, gdppos);
I_y10 = I(:, ypos);
I_pdm = I(:, pdpos);
I_divgrm = I(:, divgrpos);
I_dt = I(:, dtpos);
I_dg = I(:, dgpos);
I_coint_tax = I(:, coint_tax);
I_coint_spending = I(:, coint_spending);
I_coint_div = I(:, coint_div);

% actual, non-detrended series of tax and spending
X2t = X2;
X2t(:, coint_tax) = log(tau) - mean(log(tau));
X2t(:, coint_spending) = log(g) - mean(log(g));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate Psi and Sig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psi = zeros(N, N);
Psi_se = zeros(N, N);

for i = [inflpos ypos gdppos divgrpos pdpos dtpos dgpos]
    regr = ols(Y(:, i), [ones(T - 1, 1), X(:, :)]);
    c(i) = regr.beta(1);
    c_se(i) = regr.bstd(1);
    Psi(i, :) = regr.beta(2:end)';
    Psi_se(i, :) = regr.bstd(2:end)';
    R2(i) = regr.rsqr;
    R2bar(i) = regr.rbar;
    eps(:, i) = Y(:, i) - c(i) - X * Psi(i, :)';

    clear regr;
end

Psi(coint_tax, :) = Psi(dtpos, :);
Psi(coint_spending, :) = Psi(dgpos, :);
Psi(coint_div, :) = Psi(divgrpos, :);
Psi(coint_tax, coint_tax) = Psi(coint_tax, coint_tax) + 1;
Psi(coint_spending, coint_spending) = Psi(coint_spending, coint_spending) + 1;
Psi(coint_div, coint_div) = Psi(coint_div, coint_div) + 1;

Psi_se(coint_tax, :) = Psi_se(dtpos, :);
Psi_se(coint_spending, :) = Psi_se(dgpos, :);
Psi_se(coint_div, :) = Psi_se(divgrpos, :);

Sigma = cov(eps(:, [inflpos ypos gdppos divgrpos pdpos dtpos dgpos]));
tmp = chol(Sigma, 'lower');
Sig = zeros(N, N);
Sig([inflpos ypos gdppos divgrpos], [inflpos ypos gdppos divgrpos]) = ...
    tmp(1:4, 1:4);
Sig(pdpos, [inflpos ypos gdppos divgrpos pdpos]) = tmp(5, 1:5);
Sig(dtpos, [inflpos ypos gdppos divgrpos pdpos dtpos]) = tmp(6, 1:6);
Sig(dgpos, [inflpos ypos gdppos divgrpos pdpos dtpos dgpos]) = tmp(7, 1:7);

Sig(coint_div, :) = Sig(divgrpos, :);
Sig(coint_tax, :) = Sig(dtpos, :);
Sig(coint_spending, :) = Sig(dgpos, :);

eps = [eps(:, [inflpos ypos gdppos divgrpos]), eps(:, divgrpos), ...
                   eps(:, pdpos) eps(:, dtpos), eps(:, dtpos), ...
                   eps(:, dgpos), eps(:, dgpos)];
Sigma = Sig * Sig';
